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FOREWARD 


This  report  is  the  annual  report  of  the  research  project  entitled 
"Fluid  and  Thermodynamic  Characteristics  of  Compressible  Recoil  Mechanism" 
covering  the  period  from  June,  1978  to  August,  1979.  The  project  origi¬ 
nally  was  funded  on  February  1,  1977  for  a  three  year  period  of  research 
by  the  U.  S.  Army  Weapon  Command  of  Rhode  Island,  Illinois.  However,  the 
project  was  interrupted  after  the  first  year  on  February  1,  1978  because 
of  the  reorganization  of  the  U.  S.  Army  Armament  Research  and  Development 
Command.  The  project  was  then  reactivated  by  the  U.  S.  Army  Armament 
Research  and  Development  Command  from  Dover,  N.  J.  through  U.  S.  Amy 
Research  Offices  in  North  Carolina  under  the  contract  DAAG  29-78-G-0120 
under  a  reduced  funding  level  at  approximately  a  one-half  budget  for  the 
period  from  June,  1978,  to  August,  1979.  The  research  effort  thus  is 
accordingly  also  reduced.  The  portion  of  the  proposed  experimental 
investigation  was  thus  terminated  and  the  effort  is  concentrated  in  the 
development  of  the  numerical  prediction  for  the  fluid  motion  in  the 
recoil  mechanism.  Substantial  progress  has  been  made  in  the  numerical 
solution  to  the  Navier-Stokes  equation  at  high  Reynolds  number  from  the 
development  of  the  new  numerical  scheme  called  "Finite  Differential 
Method".  An  outline  of  the  procedure  of  applying  the  method  to  the 
recoil  mechanism  is  given  which  can  be  carried  out  in  the  third  year 
effort. 

The  project  was  conducted  by  Professors  C.  J.  Chen  and 
E.  0.  Macagno  with  graduate  research  assistants  N.  H.  Naseri-Neshat 
and  A.  M.  Fard-Behbahani. 
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I.  INTRODUCTION 

The  recoil  mechanism  is  a  component  of  the  recoil  system  which  pro- 
vides  a  retarding  force  acting  on  the  gun  so  as  to  brake  the  rearward 
motion  in  the  firing  cycle.  The  counter-recoil  mechanism  is  also  a  com¬ 
ponent  of  the  recoil  mechanisms  which  stores  the  energy  during  the  recoil 
and  dispenses  it  during  the  counter-recoil  phase  so  as  to  bring  the  gun 
back  into  its  in-battery  position.  Without  a  recoil  system,  the  force 
acting  on  the  gun  and  its  mounts  induced  by  firing  the  charge  is  so 
large  that  damages  to  the  gun  often  immediately  result.  An  improperly 
designed  recoil  mechanism  can  also  affect  the  accuracy  of  the  gun  in  hit¬ 
ting  the  target  and  reduce  the  life  of  the  gun.  Therefore,  the  recoil 
system  is  one  of  the  most  important  components  of  a  gun  system. 

Frequent  failure  and  the  undesirable  characteristic  of  many  springs 
used  in  the  recoil  mechanism  prompted  recent  designers  of  a  gun  system  to 
adopt  more  fluid-type  recoil  mechanisms  (hydropneumatic)  or  the  combination 
of  fluid  and  spring  type  recoil  mechanisms  (hydrospring) .  The  fluid  recoil 
mechanism,  gas  or  liquids,  is  throttling  during  a  firing  cycle  through  an 
orifice  which  divides  the  recoil  mechanism  into  two  or  more  chambers.  The 
difference  in  the  pressure  between  the  chambers  thus  provides  the  retarding 
force  for  braking  action.  The  retarding  force  and  the  counter-recoiling 
force  are  thus  functions  of  the  fluid  motion,  the  type  of  fluids,  and  the 
orifice  geometry  of  recoil  mechanisms.  Therefore,  the  knowledge  of 
fluid  mechanics  of  the  recoil  mechanism  is  an  indispensable  part  of  the 
weapon  design. 

The  fluid  type  of  record  mechanism  was  considered  in  the  1960's  by 
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the  Rock  Island  Arsenal  and  the  Iowa  Institute  of  Hydraulic  Research 
[1-5].  However,  it  was  found  [4]  that  the  designed  recoil  characteristic 
such  as  retarding  force  and  recoiling  motion  often  do  not  agree  with  the 
experimental  measurement.  The  reason  for  this  discrepancy  is  largely  due 
to  lack  of  detailed  knowledge  about  the  fluid  dynamic  characteristics  of 
the  recoil  mechanisms.  In  particular,  the  value  of  the  orifice  discharge 
coefficient  which  determines  the  flow  rate  through  the  orifice  and,  hence. . 
the  retarding  force  during  recoil  motion.  The  value  of  the  orifice  dis¬ 
charge  coefficient  in  many  designs  was  assumed  to  be  constants,  typically 
ranging  from  0.6  to  0.9,  which  were  obtained  from  experiments  of  steady 
flow  through  stationary  orifice.  However,  the  unsteady  flow  in  the 
recoil  mechanisms  is  created  by  the  unsteady  motion  of  the  orifice. 

Indeed,  the  orifice  discharge  coefficient  should  be  the  orifice  discharge 
function  and  should  depend  on  the  flow  velocity  and  the  location  and 


geometry  of  the  orifice. 

Several  experimental  data  for  simple  geometry  and  steady  flow  through 
stationary  orifice  are  available  [6J.  Theoretically  calculated  orifice 
discharge  coefficient  for  small  Reynolds  number  with  stationary  orifice 


are  also  available  [7].  However,  neither  experimental  data  nor  theoreti¬ 
cal  prediction  of  the  orifice  coefficient  are  available  for  the  moving 
orifice  in  unsteady  flow  with  large  Reynolds  number  such  as  in  recoil 
mechanisms. 

The  aim  of  this  project  is  to  provide  the  fundamental  description 
of  fluid  mechanics  for  recoil  mechanisms  and  to  present  an  approximate 
determination  of  orifice  coefficient  and  to  describe  the  fluid  motion 
in  the  recoil  mechanism. 
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II.  ONE  DIMENSIONAL  ANALYSIS  OF  ORIFICE 

AREA  AND  ORIFICE  DISCHARGE  COEFFICIENT 


In  order  to  examine  how  the  fluid  flow  in  the  record  mechanism  affect 
the  record  motion  and  how  the  orifice  area  is  designed,  a  one  dimensional 
analysis  of  fluid  motion  in  the  record  mechanism  ia  presented  here. 


il-l.  Criterion  of  Orifice  Design 

A  recoil  mechanism  may  be  designed  to  operate  on  incompressible  or 
compressible  principle.  For  incompressible  type,  the  total  volume  of  the 
recoil  cylinder  will  remain  the  same  during  recoil  motion  since  the  fluid 
is  treated  .-a  incompressible.  For  the  compressible  type,  the  volume  may 
be  decreased  in  the  recoil  motion.  The  fluid  used  in  the  recoil  mechanisms 
may  be  gas  or  liquid  such  as  air,  nitrogen,  hydraulic  or  silicone  oil. 

In  both  cases,  the  orifice  must  be  properly  designed  to  achieve  desired 
characteristics  of  retarding  force  for  recoil  mechanism. 

Before  one  calculates  the  fluid  motion  and  the  pressure  distribution 
in  a  recoil  mechanism,  design  characteristics  of  recoil  mechanisms  should 
be  first  defined  so  that  the  required  geometry  of  the  orifice  can  be 
determined  accordingly. 

Generally,  all  recoil  mechanisms  work  on  some  combination  of  the 
same  basic  principles;  that  of  providing  a  controlled  resistance  over  a 
set  distance  to  check  the  motion  of  the  recoiling  parts,  then  returning 
them  to  the  firing  position  and  providing  a  sufficient  restraint  to  hold 
them  in  that  position  at  a  given  elevation.  A  basic  requirement  for  a 
recoil  mechanism  is  that  the  resistance  to  recoil  should  be  nearly  con¬ 
stant  for  a  prescribed  recoil  distance,  since  this  will  produce  the 
smallest  possible  force  on  the  gun  structure.  A  typical  recoil  force 
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versus  recoil  distance  relation  is  given  in  Fig.  1.  The  area  under  the 
recoil  force-time  curve  represents  recoil  impulse  to  be  dissipated  by  the 
recoil  mechanism.  Clearly,  a  rectangular  curve  will  yield  the  lowest 
retarding  force  for  a  given  recoil  motion.  However,  a  rectangular  curve 
is  not  applicable  at  the  beginning  because  the  recoil  resistance  K(t) 
should  not  exceed  the  recoil  force  B(t)  created  by  the  propellant  gas. 

t 

The  total  resistance  force  K(t)  is  a  combination  of  a  hydraulic 
force  F0(t),  a  spring  force  if  a  spring  is  used  with  fluid  and  friction 
of  moving  mechanism.  Whichever  the  type  of  recoil  mechanism  is  used, 
this  force  works  as  a  unit. 

In  general,  the  recoil  force  B(t)  is  offset  largely  by  the  hydraulic 
resistance  FQ(t)  offered  by  throttling  the  fluid  through  the  controlling 
orifice  and  by  the  spring  force.  Thus,  design  of  the  control  orifice  is 
a  vital  part  of  the  recoil  mechanisms. 

We  now  consider  a  basic  recoil  mechanism  as  shown  in  Fig.  2  that  a 
recoiling  mass  M  is  subjected  to  a  time  dependent  recoil  force,  B(t), 
from  the  discharge  in  the  breech  of  a  gun  creating  an  acceleration  of  X 
in  the  x  direction.  Here  X  «=  0  is  taken  as  the  orifice  plane  before  the 
recoil  motion.  The  recoil  motion  of  mass  M  forces  the  fluid  to  flow  from 
Chamber  (1)  to  Chamber  (2).  The  orifice  installed  between  the  two 
chambers  restricts  the  flow  from  Chamber  (1)  to  (2)  so  as  to  create  a 
pressure  difference  between  the  two  chambers.  This  differential  pressure 
acts  on  the  wall  of  recoiling  orifice  and  A£  to  control  the  recoiling 
speed  and  motion.  The  force  balance  for  the  recoil  motion  can  be  written 
as 


Figure  2.  Typical  Recoil  Mechanism 
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MX  -  B(t)  -  K(t)  +  M  g  sin  0 

(2-1) 

K(t)  -  P2dA2  f  Pf  -  P0(t)  4-  P£ 

where  and  P2  are  respectively  the  pressure  distributed  on  the  orifice 
surface  facing  Chamber  (1)  and  Chamber  (2).  and  are  area  of  the 
orifice  facing  Chambers  (1)  a. d  (2) .  Here  A^  =  A^  is  for  incompressible 
type  and  A^  >  A2  for  compressible  type.  F^  is  tb-  frictional  resistance 
force  from  seal,  various  moving  parts  and,  if  any,  recuperator.  The 
last  term  of  Eq.  (2-1)  is  the  contribution  of  the  gravitational  force 
with  0  being  the  inclined  angle  of  the  recoil  mechanism  with  respect  to 
horizontal  plane.  The  dot  on  X  denotes  the  time  derivative.  In  Eq. 

(2-1),  the  mass  M,  the  retarding  force  from  breech  B(t),  two  areas  A^  and 
A2  and  frictional  resistance  force,  Ff,  are  assumed  known  but  distribution 
of  pressures  P^  and  P2  acting  on  the  surface  A^  and  A2  are  unknown  and  must 
be  solved  associated  with  the  fluid  motion  induced  by  the  motion  of  re¬ 
coiling  parts.  In  general,  both  P^  and  P2  must  be  integrated  from  the 
pressure  distribution  respectively  on  their  surfaces  A^  and  A^. 

Now  the  pressure  distribution  on  the  orifice  surface  and  their  time 
variations  must  be  determined  by  solving  the  fluid  motion  in  the  recoil 
fluid  chambers.  This  creates  the  coupling  between  the  recoiling  force 
and  the  fluid  motion.  It  should  also  be  pointed  out  here  that  the  de¬ 
signer  has  very  little  control  over  the  resistance  force,  F^,  from  the 

friction  and  recuperator  and  the  gravitational  force,  M  ,  but  he  can  con- 

8 

trol  the  resistance  offered  by  the  pressure  force  from  the  two  chambers 
by  properly  designing  the  orifice  such  that  the  total  resistance  to  the 
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recoil  force  is  sufficient  to  arrest  the  rearward  motion  of  the  gun  in  a 
specific  recoil  length  allowed  in  the  design. 

One  is  now  faced  with  the  question  of  what  orifice  design  will  pro¬ 
vide  the  lowest  peak  resistance  force  for  a  given  mechanism.  To  set  up 
the  criterion  for  the  orifice  design,  let  us  consider  Eq.  (2-1)  for 
recoil  motion  again.  K(t)  is  the  total  resistance  force  including  fric¬ 
tional  forces  Ff  and  pressure  forces  FQ  created  from  the  fluid  throttling 
through  the  orifice.  As  the  design  characteristic  of  recoil  mechanisms 
required  that  the  total  resistance  K(t)  to  recoil  should  be  nearly  con¬ 
stant  for  a  prescribed  recoil  length,  a  trapezoidal  shape  for  K(t)  such 
as  in  Fig,  1  is  normally  adopted. 

Since  Eq.  (2-1)  is  a  second  order  ordinary  differential  equation, 


one  needs  two  initial  conditions 

may  be  taken  as 

to  obtain  a  unique  solution. 

These  t.wo 

t  =  0 

X  =  0 

(2-2) 

t  =  0 

X* 

ii 

o 

(2-3) 

That  is  initially  the  recoil  mechanism  is  at  rest.  However,  if  the 
total  resistance  K(t)  is  not  specified,  there  will  be  infinite  many 
solutions  of  X(t)  to  Eqs.  (2-1),  (2-2)  and  (2-3).  To  determine  K(t),  we 
need  additional  design  constraints.  They  are: 

(1)  the  velocity  of  recoil  X  is  zero  at  the  end  of  the  recoil. 


£ 


1 

iii 


1 1 
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(2)  the  total  length  of  recoil  allowed  is  L,  which  must  be 


achieved  in  a  duration  of  t  ,  or 

r* 


t  -  t 


X(tr)  -  L 


(2-5) 


While  the  recoil  length  may  normally  be  specified,  the  time  duration, 
tr,  to  arrest  the  rearward  motion  of  the  recoil  is  dictated  by  the  re¬ 
sistance  force  K(t)  and  is  not  known  normally  a  priori.  Indeed,  t  must 
be  solved  with  K(t)  together.  To  determine  the  function  K(t)  we  reason 
as  follows 

First,  dC  the  beginning  of  recoil,  t  =  0,  and  at  the  end  of  recoil, 

•  *' 

t  =  tr,  the  mechanism  is  at  rest  and  zero  acceleration  X  =  X  »  0  and 
meantime  the  breech  force  is  zero  or  B(0)  =  B(ty)  =  0.  We  thus  require 


from  Eq.  (2-1) 


■  K(0)  =»  Mg  sin  9 
K2  =  K(tr)  =  Mg  sin  0 


(2-6) 


(2-7) 


Here  and  are  the  resistance  at  the  beginning  and  at  the  end  of 


recoil  motion. 


Secondly,  integrating  Eq.  (2-2)  twice  and  with  the  initial  conditions 


(2-3)  and  (2-4),  we  obtain 


f 


B(t)  dt  -  K(t)  dt  +  Mg  sin  0  t 
J0 


(2-8) 


*  f  f  B(x)dxdt  -  f  f  K(x)dxdt  +  Mg  sin  0  ~- 

J o  Jo  Jo  Jo  1 


rt  rt 


(2-9) 


Now  with  additional  design  constraints  (2-4)  and  (2-5),  Eqs.  (2-8)  and 


(2-9)  become 
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ffcr  ffcr 

0  -  B(t)dt<-  K(t)dt  +  Mg  sin  6  tfc 
■*0 


(2-10) 


and 


ffcr  fcr  f 

B(t)dtdt  - 

•  0  '0  'o  ' 


ft 


ML 


0  '0 


KCOdtdt  +  Mg  sin  0  —  (2-11) 


Integral  Eqs.  (2-10)  and  (2-11)  thus  provide  us  with  a  means  of  determining 
the  recoil  time  duration  tr  and  the  total  resistance  force  K(t).  Since 
the  total  resistance  force  K(t)  is  generally  desired  to  be  of  trape¬ 
zoidal  shape  as  shown  in  Fig.  1,  if  the  initial  and  final  periods,  t-^ 
and  t^,  of  the  variable  resistance  force  are  specified  then  Eqs.  (2-10) 
and  (2-11)  determine  uniquely  the  time  tr  and  the  maximum  total  resistance 

for  K  . 

max 

Either  the  moment  area  method  or  a  trial-and-error  procedure,  as 
discussed  by  Arora  and  Haug  [1]  and  Coberly  [2],  may  be  used  to  solve 
Eqs.  (2-10)  and  (2-11).  Readers  are  referred  to  the  above  references 
for  more  details.  The  total  resistance  function  K(t)  thus  can  be  used 
as  the  criterion  of  orifice  design. 

When  the  total  resistance  function  K(t)  is  determined  one  may  sub¬ 
stitute  the  K(t)  function  into  Eqs.  (2-8)  and  (2-9)  to  obtain  the  recoil 
velocity  X(t)  and  the  recoil  displacement  X(t)  as  a  function  of  time. 

The  recoil  velocity  and  displacement  in  turn  become  the  boundary  con¬ 
ditions  for  the  solution  of  fluid  flow  in  the  recoil  mechanism.  For  tue 
recoil  mechanism  shown  in  Fig.  2,  the  orifice  now  must  move  with 
X  velocity  in  the  recoiling  phase. 

The  fluid  throttling  through  the  orifice  must  produce  a  resistance 
force  such  that  when  combined  with  other  frictional  resistance  from 


seals  and  other  rubbing  parts,  F^,  it  is  equal  to  the  desired  total 
resistance  force,  K(t).  It  is  a  difficult  task  for  a  designer  to  design 
such  an  orifice,  since  he  must  solve  or  know  how  the  pressure  force  is 
generated  when  the  fluid  is  throttling  through  the  orifice.  What  is 
more  difficult,  the  orifice  clearance  area  in  general  is  variable  at 
different  locations;  the  flow  is  unsteady  and  may  be  in  both  laminar  and 
turbulent  regions. 

II-2.  Quasi-Steady  One-Dimensional  Solution 

The  fluid  motion  and  the  pressure  distribution  on  the  orifice 
surface  are  governed  by  the  mass  conservation  equation,  the  equation  of 
state,  the  Navier-Stokes  equation  and  energy  equation  for  unsteady 
compressible  fluid.  In  principle,  the  continuity  equation,  momentum 
equation,  equation  of  state  and  energy  equation  -  total  of  six  equations 
provide  the  solution  of  pressure,  density,  temperature  and  three  velocity 
components.  These  equations,  however,  are  coupled  with  Eq.  (2-1)  through 
the  moving  boundary  condition  of  X(t).  Thus,  the  fluid  motion  must  be 
solved  with  the  piston  moving  at  a  speed  of  X(t) . 

The  initial  condition  for  the  flow  problem  may  be  assumed  to  be  at 
rest.  Thus,  the  velocity  is  zero  everywhere.  The  initial  density  and 
pressure  in  the  recoil  mechanism  is  uniform. 

After  the  gun  is  fired,  the  boundary  condition  for  the  fluid  motion 
required  that  the  fluid  on  the  recoil  rod  assumes  the  recoil  velocity 
X(l)  and  Is  at  rest  on  the  stationary  cylinder  surface.  Now  the  orifice 
cross  section  must  be  designed  so  that  the  pressure  forces,  Fq,  created 
by  the  throttling  of  fluid  between  the  two  chambers  plus  other  resistance 
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Forces,  F^,  is  equal  to  the  required  design  characteristic  of  the  total 
resistance  force  K(t).  At  present,  no  such  general  solution  of  fluid 
motion  exists  even  for  the  simplest  recoil  mechanism. 

We  consider  a  simple  fluid  dynamic  model  for  the  recoil  mechanism. 

This  model  in  many  applications  does  not  provide  an  accurate  description 
of  the  flow  pattern  in  the  recoil  mechanism,  but  it  provides  the  first 
approximate  design  of  the  control  orifice.  This  method  also  provides  a 
basis  for  further  modifications  in  cases  of  complex  geometries  and  com¬ 
pressible  flow  and  the  flow  is  one  dimensional.  The  assumption  of  a 
quasi-steady  flow  is  equivalent  to  the  assumption  that  the  flow  is  in¬ 
stantaneously  steady  and  that  the  force  due  to  the  acceleration  and 
deceleration  of  fluid  in  the  recoil  is  negligible.  The  one  dimensional 
assumption  implies  that  the  pressure  is  uniform  in  the  radical  section  of 
the  mechanism.  The  resistance  force  offered  by  the  orifice  in  Eq.  (2-1) 
thus  may  be  written  as 

FQ(t)  -  K(t)  -  Ff  =  -  P2A2 

For  incompressible  type  of  recoil  mechanisms,  we  have  A^  =  A2  or 

(?1  -  P2)  =  Fo(t)/A1  (2-12) 

The  continuity  equations  for  one-dimensional  steady  incompressible  flow 
is  approximately  (see  Figure  2) 

Vo  “  <A1 +  VD1  (2-13) 

Hero  A  and  A,  are  respectively  the  orifice  clearance  area  and  orifice 
o  I 

surface  area  facing  Chamber  (1).  is  the  average  velocity  in  Chamber 
(1)  relative  to  orifice  velocity.  From  the  boundary  condition,  we  have 


uL  =  X(t) 


(2-14) 
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The  momentum  equation  under  quasi-steady,  one-dimensional  assumption  may 

be  integrated  from  the  Chamber  (1)  to  the  rear  end  of  the  orifice  as 

2  2 
U  P  U  P 

_i  +  li“_o+_£+AHf  (2-15) 

2  P  2  p 

Where  AH^  represents  the  loss  of  kinetic  energy  due  to  viscous  friction 
and,  if  the  flow  becomes  turbulent,  turbulent  shear.  There  are  six  un¬ 
knowns,  Pj,  P2,  PD,  U]_,  U0,  and  Aq,  in  four  Eqs.  (2-12),  (2-13),  (2-14) 
and  (2-15).  To  facilitate  the  solution  for  orifice  clearance  area  Aq, 
we  assume  that  the  pressure  P2  on  the  rear  orifice  surface  is  approximately 
equal  to  the  pressure  at  the  end  of  orifice  PQ  or 

P2  -  PD  (2-16) 

The  reason  for  adopting  this  approximation  is  as  follows.  As  shown 
in  Fig.  2,  the  one-dimensional  approximation  is  perhaps  invalid  immedi¬ 
ately  behind  the  orifice  because  from  experiments  conducted  by  Chen  et 
al.  [7]  as  shown  in  Fig.  3,  we  recognized  that  the  flow  always  separates 
at  the  rear  side  of  orifice  except  at  the  start  of  recoil  motion.  Therefore, 
for  one-dimensional  assumption  to  remain  valid  at  the  end  of  orifice,  the 
flow  should  be  considered  like  a  jet  (see  Fig.  3)  and  separates  at  the 
rear  end  corner  of  the  orifice.  Under  this  condition,  the  pressure  in  jet 
portion  of  the  flow  will  remain  approximately  PQ.  Thus,  for  no  other  means 
of  determining  ;he  pressure  on  the  rear  surface  of  orifice,  P2,  we  equal 
P2  to  B0  on  the  same  cross  section  which  is  the  closest  pressure  available 
under  one-dimensional  assumption. 

Furthermore,  the  loss  of  kinetic  energy  due  to  friction,  AHj,  is 
normally  proportional  to  the  pressure  difference  (P^  -  P2)/p  or 
AHf  -  (Px  -  P2)/p 


(2-17) 


STREAM  PATTERN 
OBSERVED 
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Thus,  combining  Eqs.  (13),  (14),  (15),  (16),  the  orifice  velocity  Uc 


can  be  written  as 


VW 2<fi - p2) /  1  -  Uti 


(2-18) 


Here  the  orifice  discharge  coefficient  CQ  is  introduced  to  account 
for  the  frictional  loss  and  also  to  serve  as  a  correction  factor  for  the 
assumptions  made,  namely,  one-dimensional  quasi-steady  flow  and  the 
assumptions  that  lead  to  Eqs.  (2-16)  and  (2-17).  More  discussions  on 
the  coefficient  CQ  are  given  in  Section  IV. 

Now  from  Eqs.  (2-13)  and  (2-14)  we  also  have 


Ao  =  A1 


AH 


(2-19) 


We  thus  can  solve  the  orifice  clearance  area  A  from  Eqs.  (2-18),  (2-19) 


and  (2-12)  as 


(P1  ~  P2)2 


+  1-1 


CD  V^2  +  1  -1 

X2  p  A. 


(2-20) 
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Equation  (2-20)  thus  provides  an  approximate  way  of  predicting  the  ori¬ 
fice  clearance  area  for  recoil  mechanisms.  It  should  be  noted  that  the 
orifice  area  AQ  is  a  function  of  orifice  surface  area  A^,  recoiling 
speed  X,  fluid  density,  p,  orifice  discharge  coefficient,  C^,  and  the 
resistant  force  FQ(t)  required  from  the  criterion  established  for  the 
orifice  design.  Since  both  FQ(t)  and  X  are  time  dependent  and  the 
orifice  moves  with  X,  the  orifice  area  Aq  must  in  general  vary  with  the 
recoil  distance  X.  The  orifice  area  kQ  can  now  be  predicted  if  the 
orifice  discharge  coefficient  is  known. 

II-3.  Orifice  Discharge  Coefficient  Cq 

The  orifice  discharge  coefficient  CQ  was  originally  introduced 
to  account  for  the  pressure  loss  in  a  steady  flow  through  an  orifice. 
Customarily,  it  is  treated  as  a  constant.  However,  in  case  of  recoil 
mechanisms  CD  certainly  is  not  a  constant  for  the  following  reasons. 
First  the  recoil  rod  starts  with  zero  speed  and  is  accelerated  to  a  high 
speed  flow.  Then  it  decelerates  to  zero  speed  and  reverses  the  direction 
of  motion  in  the  counter-recoil  motion.  In  the  process,  the  flow  starts 
from  laminar  motion  to  turbulent  flow  and  then  back  to  laminar  motion 
again.  The  Reynolds  number  based  on  upstream  gap  and  instantaneous 
average  velocity  may  vary  from  zero  to  an  order  of  10^.  The  orifice 
discharge  coefficient  certainly  should  vary  since  Reynolds  number  varies 
over  such  a  wide  range.  In  addition,  in  order  to  create  a  desirable 
characteristic  of  the  recoil  motion,  the  orifice  clearance  area  AQ  mast 
be  designed  to  vary.  In  summary,  the  orifice  coefficient  should  be 
dependent  of: 


■^**r<*i*  WH.Tfe 


t  g^ra*  o,  ■  ^*rr*gg>K 


A.  Geometry  of  orifice  and  shape  of  Chambers  (1)  and  (2) . 

B.  Time  variable. 

C.  Laminar  and  turbulent  phenomena. 

D.  The  average  velocity  across  the  orifice  and  viscosity  and 
linear  dimensions  of  orifice,  or  the  Reynolds  number. 

E.  Fluid  compressibility. 

The  precise  value  of  orifice  discharge  coefficient  for  recoil 
mechanisms  is  not  available  to  date  since  design  of  an  orifice  varys 
from  one  gun  to  another.  Even  for  the  same  weapon,  if  a  different 
charge  is  used,  the  orifice  coefficient  will  be  different  because  dif¬ 
ferent  recoil  forces  produce  different  recoiling  speeds.  Although  as  a 
rough  approximation,  constant  values  of  C^,  between  0.7  ~  0.9,  had  been 
used  [3,  4]  but  the  resulting  resistance  forces  were  not  satisfactory. 

In  addition,  the  resistance  force  generated  by  the  recoil  mechanism  was 
shown  by  Nerdahl  and  Frantz  [5]  and  Coberly  and  Frantz  [6]  to  be  sensi¬ 
tive  to  the  change  of  the  value  of  orifice  discharge  coefficient. 
Therefore,  a  variable  orifice  discharge  coefficient  should  be  used  when¬ 
ever  the  variable  value  of  is  available.  To  date  neither  experimental 
correlation  nor  theoretical  prediction  for  variable  CQ  is  yet  available. 
Chen  and  Macagno  [7]  are  currently  attempting  to  solve  the  detail  of 
fluid  motion  and  thereby  to  predict  the  value  of  from  the  numerically 
prediction  pressure  distribution.  While  the  progress  is  being  made,  we 
suggest  the  following  temporary  measure.  To  determine  the  variable 
orifice  discharge  coefficient,  we  propose  that  the  value  of  variable 
orifice  discharge  coefficient  at  any  instance  during  the  recoil  be  the 
value  for  steady  flow  with  the  flow  Reynolds  number  ,  UDA>>  corresponding 
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to  the  instantaneous  Reynolds  number  XD/v.  Here  D  is  the  upstream 
characteristic  length,  the  diameter  in  case  of  tube  flow  and  the  gap 
width  in  case  of  annular  flow.  U  is  the  average  velocity  upstream  of 
the  orifice.  This  assumption  is  consistant  with  the  quasi-steady 
assumption  made  in  deriving  Eq.  (2-21). 

Figure  4  shows  the  experimental  correlation  [8]  of  the  orifice  dis¬ 
charge  coefficient,  CQ,  versus  steady  flow  Reynolds  number  UD/v  as 
defined  by  Eq.  (2-18).  In  this  case  Aq  is  the  orifice  hole  area  ird^/4 
and  (A^  +  Aq)  is  the  upstream  area  or  irD^/4.  Although  the  geometry  of 
orifice  in  Fig.  4  differs  from  the  one  in  Fig.  2  for  recoil  mechanisms. 
Fig.  4  does  illustrate  the  qualitative  variation  of  orifice  discharge 
coefficient  with  Reynolds  number  as  well  as  the  geometry.  In  Fig.  4, 

L  is  the  thickness  of  the  orifice.  The  ratio  L/D  seems  to  have  an 


appreciable  effect  on  CD  value.  For  a  given  Reynolds  number  CD  value 
in  general  increases  with  L/U  ratio  except  at  a  large  Reynolds  number 


The  flow  over  Reynolds  number  of  1(T  is  likely  to  be  turbulent  in 
some  region  of  the  flow.  Turbulence  is  likely  to  be  present  near  both 
orifice  front  and  rear  extruding  corners  and  behind  the  orifice.  The 
flow  near  the  orifice  is  very  complicated.  The  flow  may  separate  at 
the  front  corner  and  reattach  and  then  separate  at  the  rear  corner 
again.  Even  when  the  upstream  slow  is  steady,  the  separation  phenomena 
may  still  be  unsteady.  That  is,  the  separation  bubble  may  grow  to  a 
certain  size  and  then  separate  from  the  corner  and  is  carried  downstream. 
A  new  separation  bubble  may  follow  to  form  at  the  front  corner  again 
and  the  phenomenon  repeats  to  give  a  definite  periodic  separation.  The 
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above  mentioned  phenomenon  was  observed  by  Chen  and  Macagno  [7]  in  simu¬ 
lation  of  recoil  motion  with  a  similar  configuration  as  shown  in  Pig.  3 
at  Reynolds  number  of  3752  based  on  upstream  gap  and  velocity. 

Examining  Fig.  4,  we  may  conclude  that  the  orifice  discharge 

coefficient  is  less  sensitive  tu  the  ratio  of  d/D  in  the  range  0.04 

to  0.25  than  to  the  ratio  of  L/D.  This  has  a  favorable  implication 

that  the  is  less  dependent  on  the  variable  orifice  clearance  area 

A  required  in  recoil  mechanisms  since  the  ratio  A  /(A  +  A)  is  similar 
o’  o  o 

2 

to  (d/D)  in  Fig.  4.  This  implication  will  lessen  the  complication  in 
determining  CD  for  recoil  mechanism.  The  variable  orifice  coefficient 
CD  to  be  used  in  recoil  mechanism  is  then  a  strong  function  of 
instantaneous  Reynolds  number  and  the  ratio  of  orifice  thickness  to 
the  upstream  gap  or  L/D. 

When  the  fluid  is  compressible,  the  orifice  coefficient  should 
be  modified.  This  is  discussed  also  in  the  report  by  Chen  and  Macagno 
[7]. 

In  summary,  if  the  resistance  force  offered  by  recoil  mechanism  is 
defined,  for  example,  a  trapezoidal  function  in  time,  then  a  criterion 
for  orifice  design  may  be  reduced  from  the  equation  of  motion  for  recoil 
mechanism.  The  orifice  clearance  area  Aq  may  then  be  approximately 
predicted  from  the  quasi-steady  one  dimensional  analysis  for  incompres¬ 
sible  flow.  However,  the  orifice  discharge  coefficient,  CQ,  must  be  pro¬ 
vided  from  experimental  correlation.  Since  the  orifice  area,  Aq,  in  general 
varies  with  recoil  distance,  a  variable  orifice  discharge  coefficient,  Cp, 
is  proposed  and,  as  a  temporary  measure,  is  taken  as  that  of  the  steady 
flow  with  the  instantaneous  Reynolds  number  XD/v  replacing  the  steady 
Reynolds  number. 


For  future  design  of  recoil  mechanism,  unsteady  two  dimensional 
analysis  of  flow  motion  is  proposed  so  that  detail  and  precise  pressure 
distribution  on  both  sides  of  the  orifice  surface  may  be  predicted. 

Such  an  analysis  may  be  made  quite  versatile  in  that  orifice  geometry, 
recoil  speed  and  distance  can  all  be  made  as  input  parameters. 


III.  SIMULATION  OF  FLUID  DYNAMIC  CHARACTERISTICS  IN  RECORD  MECHANISM 


In  this  section,  we  are  concerned  with  the  simulation  of  the  fluid 
dynamics  of  recoil  mechanisms  so  that  the  orifice  discharge  coefficient, 
the  recoil  motion  and  force  may  be  properly  predicted.  In  the  simulation 
experiment,  it  may  also  be  possible  to  study  the  flow  patterns,  and  other 
informations  which  are  helpful  in  understanding  the  dynamics  of  recoil 
systems.  The  possibility  of  using  different  media  and  materials  in  the 
simulation  is  examined.  The  design  configuration  is  made  for  compres¬ 
sible  fluid  recoil  mechanism,  however,  the  design  procedure  may  be 
applied  to  conventional  recoil  mechanisms.  The  immediate  goal  of  the 
design  is  to  obtain  experimentally  the  orifice  discharge  coefficient. 

In  the  design,  calculation  of  the  discharge  coefficient  is  accomplished 
through  the  measurement  of  forces  that  are  imparted  to  the  stationary 
and  moving  parts  of  the  recoil  mechanism,  and  the  measurement  of  fluid 
pressure  in  the  recoil  mechanism.  The  simulation  may  be  also  used  for 
comparison  with  the  numerical  predictions. 

III-l.  Simulation  Analysis 

The  purpose  of  a  recoil  mechanism  is  to  provide  a  retarding  force 
acting  on  the  recoiling  parts  of  the  gun  in  order  to  control  the  motion 
induced  by  firing  the  charge.  Two  types  of  recoil  mechanisms  for  large 
caliber  weapons  are  in  use:  (A)  compressible  fluid  (hydropneumatic) 
type,  and  (B)  spring-fluid  (gas  or  liquid)  type.  The  detailed  descrip¬ 
tion  of  recoil  mechanisms  are  available  in  [9]  and  [10].  In  the  present 
simulation  design,  a  single  orifice,  hydropneumatic  without  spring  is 
considered.  Figure  (5)  is  a  schematic  sketch  of  such  a  recoil  mechanism. 
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In  order  to  conduct  a  simulation  experiment,  the  similarity  para¬ 
meters  of  recoil  mechanism  must  be  simulated.  To  achieve  this,  let  us 
describe  the  similarity  parameters  by  examining  fluid  flow  and  piston 
motion.  The  fluid  and  piston  motions  are  governed  separately  by  two 
different  momentum  equations.  The  two  equations  are  not  independent. 

They  are  coupled  through  the  boundary  conditions  and  pressure  distri¬ 
bution  on  the  wall  surface.  Therefore,  in  the  simulation  analysis, 
both  fluid  and  solid  momentum  equations  must  be  examined  in  addition  to 
geometrical  similarity.  We  examine  these  similarities  in  detail. 

(A)  Geometry  Simulation:  The  geometric  simulation  can  be  carried 
out  by  taking  the  ratio  of  all  corresponding  lengths  of  prototype  and 
simulated  model  equal  to  an  arbitrary  constant.  Three  dimensionless 
geometric  parameters  are  considered  necessary  for  simulation.  They  are 
the  ratio  of  orifice  gaps  to  the  piston  flange,  6/H  (see  Fig.  5),  the 
width  of  the  orifice  to  the  piston  flange  (d/H) ,  and  the  ratio  of  piston 
flange  to  the  piston  radius  (H/R^)  or: 

»G1  ’  S/H’  BG2  ‘  m •  SG3  *  H/R2 

(B)  Fluid  Dynamic  Simulation:  The  fluid  dynamic  simulation  of  un¬ 
steady  flows  require  two  dimensionless  similarity  parameters  which  must 
be  kept  the  same  for  both  prototype  and  simulation  design.  These  two 
similarity  parameters  are  Reynolds  number  and  Strouhal  number.  The 
definition  of  Reynolds  and  Strouhal  number  is  obtained  from  the  govern¬ 
ing  equation  for  the  fluid  flowing  through  the  orifice.  The  Navier- 
Stokes  equation  for  an  unsteady,  axisymmetric,  laminar  flow  with  constant 
physical  properties  and  negligible  body  force  is: 
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(3-1) 


where 

2 

t  =  t*/T,  r  =  r*/R~,  z  =  z*/R„,  v  *  v*/V,  v  =  v*/V,  p  =  p*/p*v 
z  z  r  r  z  z 

where  V  and  R£  are  the  maximum  piston  velocity  (reference  velocity)  and 
the  cylinder  inner  radius  as  shown  in  Figure  5  (reference  length) ,  t* 
denotes  time,  z*  and  r*  axial  and  radial  coordinates,  v*  and  v*  axial 
and  radial  components  of  velocity,  and  p*,  p*  and  v  denote  density, 
pressure,  and  kinematic  viscosity,  respectively.  T  is  the  reference 
time  and  can  be  taken  as  the  period  of  the  recoil  motion.  Ng  and  Re  are 


Strouhal  and  Reynolds  numbers,  respectively  or: 

vR„  R 

Re  = -  and  N  =  — 

v  S  vx 


(3-2) 


Re  and  Ng  are  the  two  similarity  parameters  to  be  simulated  in  the 
design  analysis. 

(C)  Piston's  Motion  Simulation:  In  order  to  have  a  complete 
simulation  between  the  prototype  and  simulation  model,  it  is  obvious  the 
similar  piston's  motions  is  needed  for  both  cases.  This  is  to  say  that 
the  piston's  displacement  profile  must  be  the  same  for  the  prototype  as 
well  as  the  simulation  model.  The  governing  equation  for  the  motion  of 
the  piston  is  the  Newton's  second  law,  which  can  be  written  as: 


d  *  * 

a  ^  ~  W 


dt  ~(t) 


(3-3) 


where  and  denote  the  breech  and  retarding  forces  on  the  piston, 

respectively,  and  m  is  the  recoil  mass,  v*  is  the  piston  velocity  and 
t*  is  time.  Equation  (3-3)  may  be  made  dimensionless  with  the  following 


expressions: 
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t  *  * - .  V  =  —2-  R  * 

T  *  P  V  *  “(t) 


— (t *) 


■»  F', 


F*  * 
(t*) 


B  ’  (t)  F 

where  B  and  F  are  the  maximum  breech  force  and  the  maximum  retarding 
force.  Introducing  the  above  expressions  into  (3-3) 

dV 

*  "b  B(t>  -  Hf  F(t)  (3-4) 

The  dimensionless  parameters  N_  =  ——  and  N_  =  —•  must  be  the  same  for 

B  mV  F  mV 

the  simulation  design  to  achieve  the  desired  simulation.  To  facilitate 
the  calculation  of  Np  and  Ng,  F  will  be  calculated  from  the  following 
equation: 

F  =  ip  x  A  (3-5) 

where  Ap  is  the  pressure  difference  on  the  two  sides  of  the  piston  and 
A  is  the  piston  flange  area.  It  should  be  remarked  that  equation  (3-5) 
is  the  key  equation  for  coupling  fluid  motion  and  piston  motion.  The 
value  of  the  breech  force  B  can  be  obtained  from  the  recoil  force  that  is 
generated  by  the  propellant  used  in  a  particular  recoil  mechanism. 

The  goal  of  simulation  experiment  is  to  provide  information  about 
the  fluid  dynamic  characteristics  of  recoil  mechanisms,  in  particular 
the  orifice  discharge  coefficient.  To  determine  the  orifice  discharge 
coefficient  in  the  simulation  experiment,  one  may  measure  the  force, 

F,  that  is  on  the  recoil  system  stand  and  the  recoil  piston  motion 
during  the  experiment.  Through  equation  (3-5)  the  force  measurement 
provides  the  amount  of  pressure  exerted  on  the  walls  of  the  cylinder 
which  in  turn  will  be  used  to  determine  the  orifice  discharge  coefficient 
0^.  The  determination  of  the  orifice  discharge  coefficient  C^,  may  be 
derived  as  follows:  The  momentum  equation  (3-2)  with  one  dimensional 
approximation  may  be  integrated  along  the  stream  line  to  give  an 
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approximate  formula  as  shown  in  Eq.  (3-6)  where  AH  is  the  energy  loss 


due  to  viscous  effect. 


,2  ,  2 
vl_  +  £j.  =  :5L_  =  £_AH 
2gc  +  P  2g  p 


(3-6) 


In  equation  (3-6)  the  flow  is  assumed  to  be  quasi-steady  and  the  coordi¬ 


nate  system  is  fixed  on  the  cylinder,  p  is  the  density  of  fluid  and  g 


is  gravitational  constant,  p'  and  p  are  the  averaged  pressures  down¬ 


stream  and  upstream  of  the  orifice,  respectively,  v  and  v'  are  the 


velocities  at  the  upstream  location  and  at  the  orifice,  respectively. 


From  continuity  equation  for  the  incompressible  fluid,  we  have 


Av  =  A'v' 


(3-7) 


where  A  and  A'  are  the  piston  flange  area  and  the  orifice  area,  respec¬ 


tively.  Combining  equations  (3-6)  and  (3-7),  one  obtains  the  fluid 


velocity  at  orifice  area  as: 


v’  = 

p 


(3-8) 


where  CQ  is  introduced  to  account  for  the  energy  loss,  unsteady  flow  and 


one  dimensional  approximation.  The  variables  Ap,  v’,  CD  and  A’  may  be 


considered  to  jry  with  respect  to  time  in  the  recoil  system.  From 


equations  (3-8)  and  (3-7),  it  is  obvious  that  if  v,  v',  Ap,  A  and  A' 


are  known  CQ  may  be  calculated.  For  a  simulation  experiment  A  and  A' 


are  given.  The  fluid  velocity  at  the  orifice  v'  can  be  obtained  from 


the  measurement  of  piston  velocity  V  through  equation  (3-17) .  The 


determination  of  the  pressure  difference  Ap  may  be  achieved  by  measuring 


the  retarding  force  F  on  the  cylinder  with  equation  (3-5).  The  force  F 
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can  be  measured  such  as  by  a  strain  gauge  type  of  instrumentation  with 
an  experimental  arrangement  shown  in  Fig.  (6).  The  piston  velocity  can 
be  measured  also  during  the  experiment  by  displacement  transducer.  The 
force  F  can  be  obtained  as  follows:  Suppose  that  the  hydraulic  fluid  of 
the  simulation  model  is  evacuated  such  that  inside  the  recoil  system 
there  is  vacuum  (or  just  air).  Then  the  recoil  motion  of  the  piston 
will  be  met  by  only  the  frictional  force  from  the  seals  of  the  outer 
cylinder  and  may  be  recorded  as  curve  (a)  in  Fig.  (7).  Now,  with  the 
simulated  hydraulic  fluid  filling  the  recoil  mechanism,  the  recoil  motion 
will  be  retarded  by  the  additional  force  F  =  Ap^A.  This  force  F  may 
then  be  obtained  by  subtracting  the  frictional  force  (curve  a)  from  the 
force  received  by  the  external  cylinder  of  curve  (b)  on  Fig.  (7).  The 
difference  between  curve  (a)  and  (b)  is  due  to  the  force  exerted  by 
fluid  on  the  orifice  and  hence  the  recoil  mass.  Once  F  is  found,  Ap 
can  be  calculated  and,  using  equation  (3-8),  CD  can  be  plotted  versus 
time. 

III-2.  Simulation  Design 


As  stated  previously,  the  present  design  of  simulation  experiment  is 

for  determination  of  orifice  discharge  coefficients  in  recoil  mechanisms. 

In  order  to  simulate  the  fluid  flow  in  the  recoil  mechanism,  the  shape 

of  the  geometry  and  the  geometrical  ratios  N^  =  6/H,  N^  =  d/H,  = 

H/R  ,  Strouhal  number  Ng  and  Reynolds  number.  Re,  must  be  simulated  for 

both  the  prototype  and  the  simulation  model.  In  order  to  simulate  the 

recoil  motion,  the  parameters  Nj,  and  must  also  be  simulated  for  the 

prototype  and  the  simulation  design.  Based  on  available  date  [10]  Ng, 

R  ,  N^  and  N„  can  be  calculated  for  the  prototype.  One  of  the  simulation 
e  F  B 
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consideration  is  that  the  mass  of  the  recoil  piston  in  the  simulation 
experiment  should  be  reasonable.  Another  consideration  is  that  the 
time  duration  of  the  experiment  should  be  reasonably  long  to  permit 
accurate  measurements.  Table  1  shows  the  properties  of  water  and  Dow 
Corning  210-10  centistoke  silicone  fluid.  Table  2  shows  the  properties 
of  "clear  Incite"  which  can  be  used  for  the  cylinder  in  the  simulation. 
Clear  lucite  is  a  material  strong  enough  to  withstand  the  simulated 
stresses.  It  is  transparent,  allowing  observation  of  fluid  motion  during 
the  experiment. 

In  order  to  achieve  the  simulation,  we  choose  the  simulation  fluid 

—5  2 

to  be  water  which  has  a  smaller  kinematic  viscosity  (1.06  x  10  ft  /sec) 

than  the  hydraulic  oil(typically  10  4  ft  /sec) .  The  smaller  the  fluid 

kinematic  viscosity  is,  the  smaller  the  simulated  piston  velocity  is 

VR 

required  to  achieve  the  same  Reynolds  number _ 2.  This  is  good  because 

v 

it  leads  to  a  longer  observation  time  for  the  simulation  than  that  would 
be  in  the  real  recoil  motion. 

Consideration  of  using  water  provides  many  advantages  in  the  simu¬ 
lation.  In  the  following,  we  examine  the  use  of  water  as  the  simulation 
fluid  instead  of  the  Dow  Corning  oil.  Table  3  shows  different  parameters 
and  symbols  with  their  definition.  These  values  are  taken  from  [10]  the 
recoil  mechanism  used  in  the  simulation  analysis.  In  the  complete  simu¬ 
lation,  all  geometry  ratios  =  <5/H,  “  d/H  and  Ng^  =  H/Rj  must  be 

k’pt  similar.  We  now  need  to  simulate  geometric  similarity  Ng^,  N^,  Ng^, 
fluid  flow  similarity,  Ng,  Re  and  recoil  similarity  Np  and  Ng.  Table  4 
gives  the  characteristic  value  of  various  simulation  (see  Fig.  5)  where 


TABLE  1.  PROPERTIES  OF  FLUIDS 


Dow  Corning  Silicone  Fluid  210-10 


Specific  gravity 
Density  (at  STP) 

Coefficient  of  Thermal  Expansion 
Viscosity 

Power  Point  at  St  d  Pressure 
Flash  Point  at  St  d  Pressure 
Specific  Heat  at  77 °F 
Boiling  point  at  0.5  mm  Hg. 


Water 

Specific  gravity 
Density  (at  STP) 

Coefficient  of  Thermal  Expansion 
Viscosity 

Specific  Heat  at  77°F 


0.94  - 

0.03396  lbm/in 
0.006  in/m/°F  at  STd  P. 

10  centistoke 

-85°F 

325°F 

0.4  Bt/lbm  °F 

Absolute  Pressure  +  392°F 


0.036  lbm/in 
0.148  1/°R 
■*.42  Centipoise 
0.99  Btu/lbm  °F 


TABLE  2.  PROPERTIES  OF  CLEAR  LUCITE 


Specific  Gravity 

1.19 

Tensile  Strength 

10500  (psi) 

Flexural  Strength 

16000  (psi) 

Compressive  Strength 

18000  (psi) 

Modulus  of  Elasticity 

45000  (psi) 

Coefficient  of  Thermal  Expansions 

0.00005  (in/ in°F) 

Refractive  Index 

1.49 
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TABLE  3.  LIST  OF  SIMULATION  VALUES 


Symbol 

R2  =  Cylinder  inner  radius 

7  in 

v  =  Kinematic  viscosity 

1  07  =■  10~4  ft2/sec 

v  ■  Piston  velocity 

33.33  ft/sec 

v*  =  Fluid  velocity  at  orifice 

558.735  ft/sec 

A  =  Piston  flange  area 

57.005  in2 

A'  =  Orifice  area 

4  in2 

x  =  Period  of  a  cycle 

0.05  sec 

m  =  recoil  mass 

2855  lbm 

p  =  density 

0.03396  lbm/in^ 

gc  ■  gravity  acceleration 

32.2  lb. ft/lb.  sec2 

F  =  Retarding  force 

391271.15  lb 

B  =  Breech  force 

Ng  =  Strouhal  number  (R^/vr) 

0.35 

R  =  Reynolds  number  (R_v/v) 

“  4 

3046063 

Np  =  Fr/mv 

0.1935 

=  Piston  outer  diameter 

10.25  in 

6  =  Orifice  opening 

0.09  in 

H  -  Flange  Height 

1.775  in 

d  =  Flange  thickness 

1  in 

L  “  Cylinder  length 

18  in 

K  *  0.357  means  that  the  simulated  size  is  only  35.7%  of  the  prototype 
size. 

The  reason  for  considering  full  scale  simulation  is  that  the  proto¬ 
type  may  be  obtained  from  the  existing  recoil  mechanism  in  use  without 
the  need  of  constructing  a  new  simulating  model.  The  use  of  the  same 
hydraulic  fluid  in  the  simulation  means  that  the  simulation  is  identical 
to  the  full  scale  firing  of  a  weapon  which  is  difficult  to  perform  in 
the  laboratory.  On  the  other  hand,  if  simulation  with  actual  size  but 
with  water  as  media  is  considered,  velocity  is  one-tenth  of  the  actual 
speed  and  the  time  scale  is  10  times  larger  for  the  simulation  and  the 

simulated  recoil  force  (B)  is  reduced  ti  ibe-hundredth  of  the  real  recoil 

s 

force  (B)  .  The  simulated  mass  remains  the  same  as  that  of  the  prototype. 
P 

Therefore,  the  full  size  simulation  with  water  as  hydraulic  oil  offers 
advantages  of  (1)  slower  simulation  velocity,  (2)  longer  simulation  time 
and  (3)  smaller  simulated  recoil  force.  If  the  full  size  simulation  is 
still  undesirable  because  of  its  dimension  or  of  its  large  mass,  the 
simulation  with  reduced  size  is  calculated  and  tabulated  in  Table  4. 
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IV.  TWO  DIMENSIONAL  FLOW  IN  RECOIL  MECHANISMS  - 
A  FINITE  DIFFERENTIAL  METHOD 

IV-1.  The  Finite  Differential  (FT)  Method 

In  the  1978  annual  report  [10]  of  the  project  and  a  subsequent 
report  [11],  the  finite  difference  solution  of  unsteady  flow  in  a  recoil 
mechanism  was  given  for  a  Reynolds  number  of  3470.  Here  the  Reynolds 
number  was  defined  with  the  maximum  velocity  of  the  recoil  motion  and 
the  gap  between  the  cylinder  and  the  recoil  rod.  The  result  showed  that 
the  finite  difference  solution  with  a  grid  size  of  (20  x  50)  at  Reynolds 
number  of  5  is  stable  but  the  solution  is  unstable  at  the  Reynolds  number 
of  3470  and  needs  a  scheme  that  averages  the  vorticity  solution  to 
stabilize  the  numerical  result.  The  numerical  solution  thus  only 
qualitatively  agrees  with  the  experimental  visualization  of  the  stream¬ 
line  pattern  in  the  upper  chamber  (the  high  pressure  side)  of  the  recoil 
mechanism.  Further  survey  of  finite  difference  methods  found  that  the 
finite  difference  solution  of  the  Navier-Stokes  equation  is  quite  unstable 
at  the  high  Reynolds  number  unless  the  grid  size  is  greatly  reduced. 
Approximately,  the  grid  size  Ax  required  for  a  stable  solution  is  pro¬ 
portional  to  the  inverse  of  the  Reynolds  number.  The  increase  in  the 
number  of  grids  at  high  Reynolds  number  requires  considerably  large 
computer  storage  and  computational  time.  Even  so,  most  of  finite  dif¬ 
ference  methods  show  considerable  instability  if  the  right  combination 
of  parameters  such  as  the  time  step,  the  relaxation  factor  and  the  method 
of  solving  simultaneous  equation  were  not  adopted.  This  difficulty  also 
persists  in  many  simpler  flow  configurations  such  as  a  transient  or  a 
steady  flow  in  a  square  cavity.  For  this  reason,  a  new  numerical  method 


called  the  finite  differential  (FT)  method  was  conceived  and  developed. 

The  finite  differential  (FT)  method  is  not  the  finite  difference  (FC) 
method,  nor  the  finite  element  (FE)  method.  The  basic  principle  of  the 
finite  differential  (FT)  method  is  outlined  here  and  then  the  finite 
differential  (FT)  solution  for  the  flow  in  the  square  cavity  is  demon¬ 
strated.  A  better  convergence  and  stability  of  the  numerical  solutions 
are  obtained  from  the  finite  differential  (FT)  method. 

IV-2.  The  Principle  of  Finite  Differential  (FT)  Method 

The  basic  idea  of  the  FT  method  is  the  incorporation  of  local 
analyti '  solutions  in  the  numerical  solution  of  partial  differential 
equation  (PDE) .  The  FT  method  starts  with  subdivision  of  the  total  flow 
region,  D,  into  small  subregions  as  shown  in  Figure  8  in  which  the 
governing  partial  differential  equation  (PDE)  may  be  solved  analytically. 

The  assembly  of  all  the  local  analytic  solutions  thus  constitutes  the 
numerical  solution  of  the  problem.  Details  of  the  FT  method  are  outlined 
below. 

Consider  a  partial  differential  equation,  (PDE)  with  a  inhomogeneous 
term  G,  L((J>)  =  G,  where  L  can  be  any  partial  differential  operator, 
linear  or  nonlinear.  This  PDE  is  to  be  solved  in  a  region  D,  see  Figure 
1.  Depending  on  the  operator  L,  let  the  boundary  conditions  and/or 
initial  conditions  be  specified  so  that  the  problem  is  well  posed.  A 
numerical  solution  is  sought  when  the  problem  cannot  be  solved  analytically. 
In  order  to  solve  the  problem  with  the  finite  differential  method,  the 
complex  geometry  of  the  problem  is  broken  up  into  a  number  of  subregions 
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where  analytic  solutions  can  be  obtained.  Let  the  region  D  be  subdivided 
into  small  rectangles  by  passing  orthogonal  lines  through  the  region. 

The  intersection  of  these  lines  forms  the  nodal  points  with  I  =  1,2,3,..., 

i-1,  i,  i  +  1 .  IN,  and  J  =  1,2,3 . J-l,  j,  j  +  1,...,JN.  A  typical 

subregion  of  the  problem  with  the  node  point  P(i,  j)  may  be  surrounded  by 
the  neighboring  node  points  E  (east),  W  (west),  S  (south),  N  (north),  NE 
(northeast),  NW  (northwest),  SE  (southeast)  and  SW  (southwest),  which 
corresponds  to  points  (i  +  1,  J) ,  (i  -  1,  j) ,  (i,  j-l),  (i,  j  +  1) , 

(i  +  1,  j  +  1),  (i  -  1,  j  +  1),  (i  +  1,  j  -  1)  and  (i  -  1,  j  -  1),  res¬ 
ect  ively. 

Once  the  region  D  has  been  subdivided  into  simple  rectangular  sub- 
regions,  an  analytic  solution  in  the  single  subregion  may  be  obtained. 

In  the  case  when  the  PDE  is  nonlinear,  the  nonlinear  equation  may  be 
locally  linearized  in  the  simple  region.  However,  the  overall  nonlinear 
effect  can  still  be  preserved  by  the  assembly  of  local  analytic  solutions 
which  constitute  the  numerical  solution  of  the  PDE  over  the  whole  region 
D.  Indeed,  the  local  linearization  technique  is  also  used  both  in  the 
FC  and  the  FE  methods. 

The  problem  has  now  been  reduced  into  one  with  many  finite  regions, 
where  analytic  solutions  can  be  obtained  if  the  boundary  and  initial 
conditions  in  each  simple  finite  subregion  are  properly  specified,  thereby 
making  each  simple  subproblem  well  posed. 

Let  the  governing  equation  in  a  simple  subregion  be  L(«f> )  =  G,  where 
L  here  denotes  a  linear  partial  differential  operator  so  that  an  analytic 
solution  can  be  obtained  for  the  subregion  as  a  function  of  the  boundary 


conditions; 
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(4-1) 
(4-2) 

where  the  fN,fg,fg  and  are  prespectively  the  northern,  southern, 
eastern  and  western  boundary  conditions  of  the  subregion,  h  is  the  dis- 
tnce  between  the  point  P  and  the  sides  N  and  S.  Similarly  k  is  the  respec¬ 
tive  distances  between  point  P  and  the  sides  E  and  W.  The  boundary  con¬ 
ditions  f„  and  of  the  element  are  functions  of  x,  while  f„  and  fIT  are 

N  b  IS  W 

functions  of  y.  The  boundary  functions  f's  (i.e.,  f^,fg,fg  and  f^)  may 
be  approximately  expressed  in  terms  of  the  nodal  values  along  the  boundary 
such  as,  f's  =  1 ($„>••*)  where  <j>n,  are  the  values  of  the  dependent  vari¬ 
ables  on  each  node  points  n,  n  being  E,  W,  N,  S,  NE,  NW,  SE  and  SW  in  this 
particular  case.  Substituting  the  boundary  conditions  expressed  by 
into  Eq.  (4-2)  and  evaluating  the  relationship  between  the  functional  value 
at  an  interior  point  of  the  local  subregion  p  and  its  surrounding  points  we 
have 


L(4>)  =  G 

<f>  =  f(fN(x,t),  fg(x,t),  f£(y,t),  fy(y,t) ,  h,  k,  x,  y,  t,  G) 


<j>p  f  ^E’^W’^N’^S’^NE’^NW^SE’^SW’***^ 


(4-3) 


which  is  the  fundamental  formula  for  the  present  FT  method.  For  the 
linear  or  locally  linearized  problem,  the  9-point  FT  formula  has  the  form, 


*p  =  CE*E  +  Vw  +  CN*N  +  Vs  +  CNE*NE  +  CSE*SE 


+  CNW*NW  +  CSW*SW  +  S<G) 


(4-4) 


where  the  coefficients  C's  are  obtained  from  the  local  analytic  solution. 
S(G)  is  a  part  of  the  inhomogeneous  solution.  If  the  larger  subregion 
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shown  by  the  dotted  lines  in  Figure  1  is  used,  a  more  accurate  17-point 
formula  can  be  obtained.  It  should  be  remarked  here  that  the  finite 
differential  solution  obtained  in  Eq.  (4-4)  in  the  interior  of  the  subregion 
is  exact  in  the  sense  that  it  is  obtained  from  an  analytical  solution  to  the 
PDE  in  the  finite  subregion.  The  only  approximation  involved  is  from  the 
boundary  conditions. 

In  an  internal  finite  subregion  of  the  total  region  D,  the  neighboring 


points  of  <b  such  as  <(>_,  etc.  are,  in  general,  unknown.  However, 
p  E  W 

they  can  be  in  turn  expressed  as  an  analytic  function  of  their  neighboring 
points.  This  may  be  done  repeatedly  for  all  the  unknown  nodes  (I,J)  in 
the  total  region  D 


*i,j  ~  Ci+l,j  +  °i-l,j  *i-l,j  +  *•*  +  fi,j(G) 


(4-5) 


where  <j>.  .  is  of  a  given  subregion  and  other  4>’s  in  Eq.  (4-5)  are  the 
J  P 

boundary  values  given  in  Eq.  (4-4).  The  assembly  of  all  the  expressions 


for  all  nodes  points  can  then  be  expressed  in  a  matrix  form.  The  system 
of  algebraic  equations  can  now  be  solved  numerically  as  in  the  finite 
difference  method  to  give  the  numerical  solution  of  the  total  problem. 
This  is  the  essence  of  the  finite  differential  (FT)  method. 

There  is  an  essential  difference  between  the  finite  differential 
(FT)  method  just  described  and  the  finite  difference  (FC)  or  finite 
element  (FE)  methods.  In  the  FC  method,  the  relationship  of  4> to  its 
neighboring  points  is  not  obtained  from  the  analytic  solution  of  the 
differential  equation  but,  instead,  from  the  difference  formula  truncated 


from  the  Taylor  series  expansion  of  the  dependent  variable  about  its 
neighboring  points.  On  the  other  hand,  the  FE  method  assumes  an  approx¬ 
imated  functional  form  (shape  function) ,  normally  some  polynomial  of  lower 
degree,  say  up  to  5th  or  6th  degree,  to  represent  the  solution  and  uses 


the  variational  or  Galerkin  type  of  integration  on  the  differential 
equation  to  find  the  relation  between  (p^  and  its  neighboring  points.  It 
should  be  remarked  here  that  the  derivatives  of  the  solution  from  the 
finite  differential  solution  may  be  obtained  by  differentiating  Eq.  (4-1). 
The  resulting  equation  will  provide  the  derivative  at  any  given  point  in 
the  subregion  without  loss  of  accuracy  due  to  differentiation. 

IV-3.  The  Finite  Differential  Solution  of  Unsteady 
Two-Dimensional  Namier-Stoke  Equation 

The  unsteady  two-dimensional  incompressible  viscous  flow  can  be 
formulated  with  the  stream  function  ip  and  the  vorticity  5,  as  a  dependent 
variable.  The  independent  variables  are  two  dimensionless  space  coordi¬ 
nates  (x,y)  and  the  dimensionless  time  t.  These  independent  variables 
are  normalized  with  a  reference  length  L  and  a  reference  time  scale  L/U, 
where  U  is  the  reference  velocity.  The  Reynolds  number  thus  is  defined  as 
UL/v.  v  is  the  kinematic  viscosity.  In  order  that  the  stream  function 
satisfies  the  continuity  equation,  we  have  the  velocity  components  in  an 
x  and  y  direction,  u  and  v,  as 

u  =  \p  v  =  -ip  .  (4-6) 

y  x 

The  scalar  component  of  vorticity  (z  direction)  £  can  be  defined  as 

£  =  V  -  u  =  -(ip  +  ip  )  .  (4-7) 

*  x  y  xx  Yyy 

Taking  the  curl  operation  of  the  momentum  equation,  we  have  the  voriticy 
transport  equation 


«t  + 


+  2BE 


5  +  K 

xx  yy 


(4-8) 
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with 

A  =  R  u/2  and  R„  v/2  . 

e  e 

With  the  appropriate  boundary  and  initial  condition  equations  (4-6) ,  (4-7) > 
and  (4-8)  will  provide  the  solutions  to  u,  v,  ^  and  5*  The  pressure 
distribution  P  can  be  obtained  from  the  integration  of  the  Navier-Stokes 
equation. 

Since  equation  (4-8)  is  a  nonlinear  equation,  the  analytic  solution 
to  the  whole  region  of  flow,  either  in  a  cavity  or  in  a  recoil  mechanism 


is  not  available.  Therefore,  the  finite  differential  (FT)  numerical 

solution  is  one  way  to  obtain  an  approximate  solution.  To  implement  the 

FT  method,  the  flow  region  D  (see  Fig.  4)  is  subdivided  into  many  sub- 

regions.  In  each  region,  equation  (4-8)  is  linearized;  that  is,  A  and  B 

are  assumed  as  known  averaged  values  within  the  subregion.  In  each 

subregion,  the  boundary  conditions  at  the  north,  south,  east  and  west 

are  denoted  as  ,ip  ,ib  ,  and  €»€>€»€•  If  the  initial  condition 

n  s  e  w  n  s  e  w 

for  the  vorticity  5  is  given,  the  problem  in  the  subregion  is  well  defined 
and  can  be  solved. 

Details  of  the  finite  differential  solution  to  the  vorticity  trans¬ 
port  equation  (4-8)  and  Poisson  equation  (4-7)  are  given  respectively  in 
Appendices  A  and  B. 

The  finite  differential  solution  for  the  Poisson  equation  (4-7)  is 
given  in  Appendix  B,  equation  (B-16)  as 
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P 


00 


1 

n=l 


(c6  +  c2  +  c)  Sinf  +  I  c4 

m=l 


Sin 


mu 

2 


(4-9) 


Figure  8.  Subre 


The  velocity  components  u  and  v  may  be  obtained  by  differentiation  of 
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the  stream  function  ty.  Nine-point  formulas  for  and  are  given  by 
equations  (B-17)  and  (B-18)  of  Appendix  B. 

The  coefficients  c2,c4,Cg,  and  c  are  given  in  Appendix  B.  They  con¬ 
tain  the  known  vorticity  in  the  subregion  and  the  8-nodal  point  boundary 

stream  functions  ,\(i  tip  $  and  ).  Equation  (4-9)  is 

n6  g  sc  n  uv  w  sw  s 

thus  reduced  to  the  form  given  in  equation  (4-4)  or,  more  generally, 
equation  (4-5)  relating  the  central  nodal  value  to  the  neighboring  8 
nodal  values. 

The  steady  flow  solution  of  the  vorticity  transport  equation  (4-8) 

is  solved  to  provide  the  FT  9-point  formula  for  the  vorticity  at  the  center 
s  s 

of  the  subregion  =  £^(0, 0)  (see  Appendix  A,  equation  (A-J?),» 


Sp  -  I  <cx  +  c2)  Sin  f-  +  1  (c3  +  c4)  Sin  ^ 


(4-10) 


Again,  equation  (4-10)  may  be  written  in  the  form  given  by  equation  (4-4) 
and  (4-5). 


IV-4.  Numerical  Calculation  and  Discussion 


In  order  to  examine  the  accuracy  and  stability  of  the  new  Finite 
Differential  Method,  the  case  of  steady  recoil  motion  without  orifice  is 
considered.  This  is  the  case  similar  to  the  steady  state  driven  cavity 
flow.  The  reason  for  selecting  this  flow  is  that  there  is  an  abundance 
of  similar  works  available  for  comparison. 

To  solve  the  steady  flow  without  orifice  as  shown  in  Figure  9, 
equations  (4-6) ,  (4-7)  and  (4-8)  must  be  solved  with  the  boundary  conditions 


Recoil  Mechanism 


Figure  9  Cavity  Flow  and  Recoil  Mechanism 


•PS*! 
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x  ■  0  and  x  =  1  ,  u*v*0ori*»*0  ,  9i p/ 9x  ■  0 

y  “  0  ,  u  *  1,  v  *  0  or  f »  0  diji/dy  *  1  (4-11) 

y  *  1  ,  u  *  v  ■  0  or  tp  *  0  ,  9ip/9y  »  0 

The  finite  differential  solution  of  the  problem  is  solved  with  the  finite 
differential  9-point  formulas  (4-9)  and  (4-10)  for  the  stream  function  and 
vorticity  function.  Equations  (B-17)  and  (B-18)  in  Appendix  B  provide  the 
finite  differential  solutions  to  the  velocity  components  u  and  v.  The 
numerical  solution  procedure  is  as  follows: 

(1)  A  guessed  vorticity  is  made  and  the  stream  function  is  solved  from 
equation  (4-9)  with  the  boundary  condition  ip  *  0  given  in  equation 
(4-11) . 

(2)  The  velocity  components  u  and  v  at  each  node  are  calculated  from 
equation  (4-10)  with  the  boundary  condition  (4-11) . 

(3)  The  vorticity  function  is  then  solved  from  equation  (4-10)  with  the 
boundary  condition  (4-11) . 

(4)  The  newly  computed  vorticity  function  is  compared  with  the  guessed 
or  old  vorticity  given  in  Step  1.  The  solution  converges  if  the 
difference  is  within  the  given  convergence  criterion.  If  not,  the 
procedure  is  repeated  from  Step  (1)  using  the  newly  calculated 
vorticity  as  the  guessed  value. 

The  above  procedure  is  essentially  an  iterative  scheme.  In  solving 
each  equation,  (4-9)  and  (4-10)  an  ADI  (Alterative  Direction  Implicit) 
method  is  used  with  the  SOR  (Succession  Over  Relaxation)  scheme.  Also, 
before  returning  to  Step  1,  the  SUR  scheme  is  used  to  reduce  the  sensi¬ 
tivity  of  the  system  of  the  algebraic  equations.  Before  calculating  the 
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cavity  flow  problem  a  test  was  made  to  solve  a  simple  problem  with  constant 
velocity  u  and  v  (i.e.,  assuming  A  and  B  in  equation  (4-8)  to  be  constant). 
It  was  found  that  the  FT  solution  is  more  accurate  than  the  finite  differ¬ 
ence  method.  In  particular,  the  FT  method  showed  to  be  more  stable  than 
the  finite  difference  method  at  high  Reynolds  number  (i.e.,  large  value  of 
A  and  B) .  For  the  solution  at  high  Reynolds  number,  the  FT  method  provides 
converged  solutions  when  the  grid  size  is  reduced.  The  Reynolds  number  of 
100,  400,  and  1,000  have  been  calculated.  As  an  example,  the  cavity  flow 
with  Re  ■  1000  is  given  in  Figures  10  and  11.  Figure  10  provides  the 
stream  pattern  of  the  steady  flow  in  the  recoil  mechanism  without  orifice. 
Figure  11  gives  the  corresponding  vorticity  distribution  in  the  cavity. 

IV-5.  Conclusion 

In  this  chapter,  the  new  numerical  method  entitled  "The  Finite 
Differencial  Method"  is  introduced.  The  finite  differential  (FT)  method 
is  shown  to  give  more  accurate  results  than  the  finite  difference  (FC) 
method  under  the  same  flow  conditions.  Specifically,  the  convergence 
and  the  stability  of  the  FT  method  is  much  superior  than  the  FC  method. 

For  the  flow  with  high  Reynolds  number  of  the  order  10-3  the  FT 
method  may  still  predict  stable  solution  provided  that  small  grid  sizes 
are  used  and  the  successive  under  relaxation  (SUR)  scheme  is  used  between 
interactions. 


Although  the  FT  method  requires  much  analytic  work  in  the  initial 


I 

I 


s 

derivation  of  the  finite  differential  formula,  the  many  flow  calculation  | 

3: 

2r 

can  be  solved  relatively  easily  afterward  since  the  finite  differential  I 


formula  are  similar  even  for  different  flow  conditions.  The  FT  method 


is  ;-<)SZ  developed  and  requires  further  study  before  it  may  put  into 
practical  prediction  of  the  orifice  coefficient  c d  for  design  of  recoil 
mechanism. 


V.  SUMMARY,  CONCLUSION  AND  SUGGESTIONS 


The  present  art  of  the  design  of  recoil  mechanisms  still  relies 
heavily  on  the  one  dimensional  analysis  of  flow  motion  in  the  recoil 
system  and  on  the  trial  and  error  method  of  experiment.  The  present 
report  attempts  to  improve  the  one  dimensional  fluid  flow  analysis  in 
Chapter  II  where  the  unsteadiness  of  the  flow  in  record  mechanism  is 
approximately  incorporated  into  the  use  of  the  orifice  discharge 
coefficient.  On  the  other  hand,  a  simulation  design  is  given  in 
Chapter  III  where  an  analyses  is  given  to  show  how  the  experimental 
duration  may  be  increased  based  on  similarity  argument,  thereby  providing 
a  longer  time  duration  for  the  instrumental  measurement  as  well  as  the 
flow  visualization.  It  should  be  pointed  out  that  the  physical  under¬ 
standing  of  the  fluid  motion  in  the  recoil  mechanism  is  still  lacking. 

For  example,  design  engineers  still  do  not  know  when  and  where  the  flow 
in  the  record  mechanism  will  become  turbulent. 

In  Chapter  IV,  a  new  numerical  method  is  developed  with  the  aim  to 
calculate  fluid  motion  at  high  Reynolds  number  in  tne  record  system. 

The  preliminary  calculation  shows  that  the  FT  method  is  superior  to  the 
FC  method  in  the  accuracy  of  the  convergence  and  the  convergence  and  the 
stability  of  the  numerical  prediction.  It  is  hoped  that  the  design  of 
recoil  mechanism  may  eventually  be  done  first  on  the  computer.  Simulation 
before  an  expensive  experiment  and  the  long  duration  of  the  experimen¬ 
tation  commences.  The  modernization  of  the  design  of  the  fluid  type 
recoil  mechanism  needs  the  accurate  solution  of  the  fluid  flow  inside 
the  recoil  mechanism.  In  order  to  achieve  this  and  it  requires  a 
commitment  of  the  Array  Armament  Research  and  Development  Command. 
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APPENDIX  A.  THE  9-POINT  FINITE  DIFFERENTIAL  SOLUTION 
OF  THE  VORTICITY  TRANSPORT  EQUATION 

With  A  =  4  R  u  and  B  =  4  R  v  assumed  to  be  known  constants  in  a 
2  e  2  e 

given  subregion,  the  linearized  vorticity  transport  equation  for  a  two 


dimensional  incompressible  viscous  flow  is 


£  +  2A  ?  +  2B  K  -  5  +  S 

t  x  y  xx  yy 


(A-l) 


The  initial  condition  or  any  condition  for  the  next  time  incremental 


step  in  a  subregion  may  be  approximately  taken  as 


=  =  =  -  =2=2  =  2 

£(x,y,t  =  a  +  ay  +  a„x  +  a_xy  +  a.y  +  acx  +  a,xy 

0  12  0  4  j  b 

,  =  2  ,  =  2  2 
+  a?x  y  +  a8x  y 


(A-2) 


where  9  coefficients  a^(i  =  1,1, 2,..., 8)  can  be  determined  from  the  9 
nodal  values  of  £  at  the  boundary. 

The  boundary  condition  of  the  subregion  may  be  approximately  taken 
as,  for  example. 


SE(y>t)  =  [(aE  +  aEt)  +  (b£  +  b£t)y  +  (c£  +  cEt)y  ]  (A-3) 


Similarly  for  ?g,  and 

Introducing  a  change  of  variable 


5  =  5  e^ 


we  may  reduce  Eq.  (A-l)  to  a  simpler  form 


+  (ft2  +  b2)«  *  «*x  +  S 


(A-5) 
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with  the  initial  and  boundary  conditions 


I.C. 

t  =  t 

0 

5  »  «<*,y,t0>  e-(A*+By) 

B.C. 

x  =  Ax 

S  -  5E(y,t)  e-(Mx+8^ 

x  -  -Ax 

5  -  {w(y,t)  e<Mx-By) 

(A-6) 

y  =  Ay 

5  -  S„(x,t)  e-(Ax+B4y) 

y  =  -Ay 

5  =  «s(x,t)  a(-Ax+B4y> 

The  problem  of  (A-2)  with  (A— 6)  may  be  solved  by  superposition  of  a 

**S  — £ 

steady  part  £  and  an  unsteady  part  £  ,  or 

l  =  £S'-,y)  +  ^(x.y.t)  (A- 7) 

or 

5  =  £S(x,y)  +  Ct(x,y,t) 


For  the  FT  9-point  formula  which  provides  an  algebraic  relation  between  the 
central  node  (x  =  0,  y  =  0)  value  with  the  neighboring  8-point  node 
value  Eq.  (A-7)  can  be  evaluated  at  x  -  0,  y  =  0.  Details  of  the  derivation 
of  the  9-point  formula  is  given  below. 


(A-8) 


with  boundary  conditions 
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Is  =  5|(y)  e-(MxHiy)  ,  x  »  4x 


Is  -  5yCy)  eA4x-By  , 


x  *  -Ax 


(A-9) 


?s  rs,  ,  -(ax+BAy)  * 

K  =  5n(x)  e  3  ,  y  =  Ay 


fs  =  C|(x)  e'As+Biy  ,  y  -  -Ay 


This  problem  can  be  solved  analytically  by  further  dividing  it  into 
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simpler  problems  having  one  or  two  nonhomogenous  boundary  conditions  in 
one  direction.  Due  to  linearity,  the  final  solution  will  be  the  super¬ 
position  of  all  these  simple  problems.  Assuming  a  second  order  polynomial 
(higher  order  polynomials  can  also  be  taken)  for  boundaries  5® 

g 

and  £s>  the  coefficients  can  be  found  in  terms  of  the  surrounding  vortici- 
ties.  For  example,  for  the  east  boundary  for  the  9-point  finite  differ¬ 
ential  (steady  part  of  Eq.  (A-3)) 


s  H 
1  1 


?g(y)  =  a£  +  b£y  +  c£y2 


(A-10) 


'3  $ 

3  % 


where 


aE  5i+l,j 


vs  #• 

a  i 


bE  “  2Ay  (5i+l,j+l  "  5i+l,j-l* 


°E  "  2 Ay 2  (?i+1»J+1  "  25i+l.J  +  ^i+l.j-l5 
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Note  should  be  taken  that  if  only  a  steady  problem  is  being  solved,  the 
vorticities  in  the  above  coefficients  are  the  same  vorticities  as  being 
calculated.  However,  for  the  case  of  unsteady  Navier-Stokes  equation,  the 
vorticities  in  these  coefficients  a„,  b_  and  c0  are  the  previous  time  step 
vorticities.  In  this  case  the  unsteady  solution  is  only  part  of  the  overall 
solution  which  is  added  to  the  unsteady  part.  If  the  problem  is  calculated 
from  fluid  at  rest,  the  contribution  of  the  steady  part  for  the  first  time 
step  will  be  zero. 

The  steady  vorticity  transport  Eq.  (A-8)  is  now  divided  into  two 
simpler  problems  each  with  two  nonhomogenous  boundary  conditions  as 
follow: 

Problem  (1) 


*  0 


i.C.  L  =  (a_  +  b  v  +  c_y2)  e“(MxfBy),  x 
1  E  E*  E 


Ax 


(A-ll) 


Problem  (2) 


h  -  <aw  +  V  +  V2)  e“X'By 


x  *»  -Ax 


(A-12) 


-  0  , 


y  =  ±Ay 


u2=  0 


(A-13) 


B.C. 


l2  =  0  , 


x  =  ±Ax 


Z2  =  (^  +  bNx  +  Cj^x2)  e~(kx*&Ly'>  t  y  =  Ay  (A-14) 

^2  =  ^as  +  bsx  +  csx^  e  Ax+BAy  >  y  =  -Ay 
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where 


V 


and 


L  = 


(A2  +  B2) 


Solution  to  Problem  (1) 

Assuming  =  X(x)  Y  (y) ,  and  substituting  in  the  differential  equation, 

the  variables  are  separated.  The  two  boundary  conditions  in  the  y-direction 

in  this  case  are  used  to  find  the  eigenvalues.  Then  the  problem  is  reduced 

to  find  the  function  A  (x)  in  the  series  solution 

n 

00 

5®(x,y)  =  l  Aq(x)  sin  ~  (y  +  Ay)  (A-15) 

n=l  y 


To  do  this,  the  above  equation  is  differentiated  with  respect  to  x  and  y 

and  substituted  in  the  governing  differential  equation  (A-ll)  to  produce 

an  ordinary  differential  equation  for  A  (x)  or 

n 

A^(x)  -  (A2  +  B2  +  X2)Aq(x)  -  0  (A-16) 


When  is  the  eigenvalue  nx/2Ay.  The  two  nonhomogenous  boundary  conditions 
in  the  x-direction  are  now  used  to  find  the  solution  for  this  equation. 

The  solution  to  Problem  (1)  becomes 
00 

I^(x,y)  =  l  (C;LeEx  +  c2e  Ex)  sin  (y  +  Ay) 
n=l  y 


(A-17) 
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where 

A  (Ax)e^X  -  A  (-Ax)e 

p  —  ^  ^ 

1  “  2E  x  -2E  x  . 

e  -  e 

A  (-Ax)e^X  -  A  (Ax)e  ^X 
r  _  _n _ n _ 

2  ~  2E  x  -2E  x 

e  -  e 

with 

-AAx 

An(Al°  -  V"  (3EE1  +  bEE2  +  ceV 

AAx 

A„<-4x>  ■  v  (Vi  +  Va  +  ■=*> 

E  "  ^A2  +  B2  +  A2 


(A-19) 


(A-19) 


Eqs.  (A-13)  and  (A-1A) 
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y  *  Ay 


V  -  5„(x,t)  .-(****’> 


y  =  -Ay  ■  Sg(x,t)  e  Ax+BAy 


(A-25) 


where 


l,j..jL.jL+  {a2  +  b2, 

St  9x2  3y2 


In  order  to  solve  the  time  dependent  problem  (Ax-23)  with  given  boundary  and 
initial  conditions  (A- 24)  and  (A-25) ,  it  is  divided  into  five  simpler 
problems  each  having  one  nonzero  boundary  or  initial  condition.  These  are 
given  as  fc”„ows: 


Problem  (3) : 


LS0-° 


i.c.  cQ  -  t”_1(x,y)  , 


t  -  t0  -  o 


B.C.  5q  =  0  f 


x  *  ±Ax  (A-26) 


?0  =  °  • 


y  =  ±Ay 


Problem  (4) : 


Ui  =  0  , 


I.C.  s  =  0  , 


t  -  0 


B.C.  ^  =  (a£t  +  b£ty  +  Cgty2)  e~(AAxf0y)  f 


x  =  Ax  (A-27) 


-  0  , 


x  =  -  x 


^  °  0  , 


y  =  y 


LC2  '  0 


I.C.  =  0  * 


t  -  0 


B.C.  5,  *  0  , 


x  ■  Ax 


\  =  ^  +  bWty  +  eAAx-By  ,  X  *  -Ax 


£, -  0  » 


y  =  ±Ay 


Problem  (6] 


L?3  =  0 


«3"°  ’ 


t  =  0 


z3  *  o  , 


x  -  ±Ax 


S  =  <V  +  bNty  +  vy2)  eMx_By  , 


y  =Ay 


s3  -  0  , 


y  »  -Ay 


Ptoblem  (7) 


«4  -  0 


S4  =  0  . 


t  =  0 


5a  -  0  • 


x  =  ±Ax 


5a  =  0  , 


y  =  Ay 


C4  =  (agt  +  bgtx  +  cstx2)  e"As+BAy  , 


y  *  -Ay 
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In  the  above  problems,  the  coefficients  a,  b,  and  c  are  given  in  terms  of 
the  boundary  vorticity.  For  example,  for  the  east  boundary; 


a 

b 


c 


E 

E 

E 


5i+l,j/At 


(5 


i+lj+l 


a 


i+i,j+i 


4+i,j-i)/2AtAy 
24+l,j  +  5i+l,j-i^?AtAy 


(A-31) 


Solution  to  Problem  (3) 


Because  of  the  linearization  (assumption  of  constant  A  and  B  in  a 
subregion),  the  variables  can  be  separated  and  solved  accordingly.  The 
vorticity  £n  ^(x,y)  in  the  initial  condition  is  just  the  value  of  vorticity 
■fn  the  subregion  at  initial  time  (previous  time  step  is  denoted  by  a  super¬ 
script  n-1)  which  can  be  approximated  by  a  polynomial  of  x  and  y.  For  a 
9-point  finite  differential,  the  polynomial  and  the  coefficients  are: 

K  (x.y)  =  aQ  +  a^y  +  a£x  +  a3xy  +  a^y  +  a5x  +  agxy  +  a?x  y  +  agx  y 


where 


=  n-1 

ao  5i,j 


(A-32) 


—  (l 

.n-1 

n-i  \ 

4,j-i; 

2Ay  V 

4,  j+1 

-1-  (i 

.n-1 

n-i  \ 

2  Ax  V 

4+1,  j  " 

4-1,  j/ 

1 

/-n-1 

.  -n- 

4AxAy 

\i+l, j+1  4' 

n-1  n-1  \ 

4-1, j+1  4+1,  j-1/ 
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“S  (Cl4-1  +  Cl  1  +  25i’J) 

■5A«Z  \  i>J-L  i.J/ 


2Ay 


— — —  ( E11-^  4-  fn~^  4-  2£n~^ 

2Ax2  V^i+i.j  +  4-i,j  +  2h,i) 


a6  = 


4AxAy 


a?  =  /a  2a 
4Ax  Ay 


/-n-l  n-1  n-1  n-1  ,.n-l  or11"1  \ 

\i+l»j+l  i-l»  j**l  “  H+1J-1  "  ZH+l,j  4Vl,j/ 

/r11-^  4-  rn~^  -  £n~^  -  cn~^  _  2  rn"'^  4.  2£n*”^  ^ 

V4+i,j+]  H-1,3+1  4i+i,j-i  ”i+i,j-i  h,j+i +  Ai,j-i/ 


,  _n-l  .  _n-l 

+  5j.  1  ixi  ■»•  5 


j+1  -i-lj+l  '■i+l.j-l  i-l,j-l  *i,j  si,j+l 


/n-i 

\i+l» 

-  2£n~^  -  2En~^  -  2£n~^  ^ 

^1A-1  h-1,3/ 


+  -  2  5?'1 


(A-33) 


The  analytic  solution  to  Problem  (3)  is  found  to  be 


00  CO 


«0  -  I  I  {•'■*  - 1 +  S7  [» -  +  «2«2>  -  f  C*V5 

k=l  £=1  v 


-Tit 


„  j  kir  £ir 

x  sin  y  sin  y 


} 

(A-34) 


where 


"  =  a2  +  e2+  (Is)  +(w)3 


H  =  a.E'E'  +  a.ElE'  +  a^E'EI  +  a^ElEi  +  a.ElE'  +  aJS'El  +  a^ElEl 
Oil  121  212  322  431  313  b  3  2 


+  i?E*E*  +  SgEjE* 


Ei  = 


Ay 

e“By  sin  (y  +  Ay)dy 
J-Ay  y 


S 


Ay 

-Ay 


ye_By  sin  ~  (y  +  Ay)dy 
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■sr^ 


|>Ay 

•^-Ay 


y2  e'By  sin  —  (y  +  Ay)dy 


(A-35) 


E|,  E^  and  are  similar  to  E^,  E^  and  respectively  with  y,  Ay,  B  and 
£  changed  to  x,  Ax,  A  and  k. 

c^,,  c^,  c^  and  are  the  same  as  the  steady  part.  and  are 


Ql  =  e 


-EAx 


kir  ,,  /  i\k  2EAX-, 

11  -  e  1 

“2+fe)2 


Again  and  are  similar  to  and  with  Ax  and  k  changed  to  Ay  and 
£  respectively. 

Solution  to  Problem  (4) 

This  problem  is  solved  similarly  by  the  method  of  separation  of 

/v 

variables.  Assuming  =  F(t,x)  Y(y),  and  substituting  into  the  differential 
equation  (A-27),  the  variables  are  separated.  The  two  zero  boundary  con¬ 
ditions  at  y  =  ±Ay  are  used  to  find  the  eigenvalues.  Thus  the  problem  is 
reduced  to  finding  F^(t,x)  in 


To  find  F^(t,x),  first  a  change  of  variable  is  made  so  that  a  homogenous 
boundary  condition  in  the  x-direction  is  found  by  assuming 


h’h-  V  <SEEi  +  Vj  +  V3U 


(A-37) 


by  substitution  of  Eq.  (A-37)  into  Eq.  (A-27)  a  new  nonhomogenous  differ¬ 
ential  equation  with  the  variables  t  and  x  will  be  found  with  homogenous 
boundary  and  initial  condition  which  can  be  solved  again  by  separation  of 
variables  with  the  result 


\  ■  I  I  \  <t)  sln  ^  +  Ax>  sin  jh  (y  +  Ay)  +  l  ~2+AxA?* 

11=1  k=l  1  y  &=1 


+  sin  (y  +  Ay) 


(A-38) 


which  can  be  evaluated  at  the  center  of  the  subregion  x  =  0  and  y  =  0. 
The  coefficient  A^  (t)  is 


\(t)  "n  (t-n+£T!)T  +  ^a-e",t)6 

-*2-2  +  (t)Z  +  (^)2 

^-^?[a2  +  “2  +  fe)1  m 


,.u  ‘-•tv:. 
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Pkl  -  -24*(-l)k/  (|£) 

-AAx 

« = V~  (vi +  vj +  yp  (a-39) 

5  sIr 

aE  At  4+1, j 

^E  =  2AtAy  '5i+l,j+l  "  ^i+l.j-l5 

CE  =  ^^2  (?i+l,j+l  "  25i+l,j  +  5i+l,j-l) 

Solution  to  Problem  (5) 

The  solution  procedure  for  this  problem  is  similar  to  Problem  (4) 
but  only  the  nonhomogen  ity  is  changed  from  the  position  x  =  Ax  to  x  =  -Ax. 
The  solution  is 


00  CO 


c2  -  I  I  A,  (t)  sin  “  (x  +  Ax)  sin  (y  +  Ay) 
Jt=l  k=l  2  Z  A y 


lit 


r  x-Ax.  ^  Jin,  ,.s 

-  /  — ~ -  at  siit  -rr—  (y  +  Ay) 

L  2 Ax  2Ay  1 


(A-40) 


where 


\  (t) 
2 


1  /.  1  .  e  11  \  ,  1  ,,  -nti 

=  -  It - + - J  y  +  -  (1  -  e  )< 

n  \  n  n  /  n 


* 


n 


i 


6  =  ——  Pk2 
2Ax 

_  4Ax2 
Pk2  =  ktr 


a 


AAx 

V  + Vi + v? 


£ 


67  i 

8 


The  solution  to  Problems  (6)  and  (7)  are  similar  to  Problems  (4)  and  (5) 
with  x,  Ax,  y,  Ay,  k,  A,  A  and  B  changing  to  y.  Ay,  x.  Ax,  A,  k,  B  and  A 
respectively.  Also  the  boundaries  change  respectively  from  East  and  West 
to  North  and  South. 
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APPENDIX  B.  THE  9-POINT  FINITE  DIFFERENTIAL  SOLUTION 
OF  THE  TWO-DIMENSIONAL  POISSON  EQUATION 


The  two-dimensional  Poisson  equation  and  its  boundary  condition  here 


are  given  as 


4 + 

3xz  3y 


with  boundary  conditions 


(B-l) 


♦  *  ^E(y)  .  x  -  ax 


*  ■  *w(y)  .  x  =  -ax 


♦  =  ,  y  -  Ay 


♦  "  <l»c(x)  »  y  =  -Ay 


(B-2) 


where  for  the  9-point  finite-differential  method  the  boundaries  will  be 
second  order  polynomials  of  x  or  y.  As  an  example 


where 


*E(y)  =  aE  +  bEy  +  CEy 


aE  =  *i+l,j 


(B-3) 


bE  =  2 Ay  **i+l,  j+1  “  ^i+1, j-1^ 


2  (*i+l,j+l  "  ^i+l,  j  +  *i+l, j-l} 


To  solve  the  problem  (B— 1)  and  (B-2) ,  it  is  divided  into  two  simpler 
problems  as 
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where 

Problem  (1) : 


v2^  =  0 


II 

W 

X  = 

Ax 

N-X 

* 

11 

X  = 

-Ax 

II 

2! 

X 

y  * 

Ay 

-  i|»s(x)  , 

y  = 

-Ay 

Problem  (2) : 

“  ~S(x;y) 

<1*2  =  0  >  x  =  ±Ax 

^2  =  0,  y  =  ±Ay 


(B-4) 


(B-5) 


(B-6) 


B-l.  Solution  to  Problem  (1) 

Again  for  simplicity  and  due  to  linearity,  this  problem  is  divided 
into  two  problems  each  having  two  homogeneous  boundary  conditions  as: 


3 


if 

If 


^(x.y)  =  u(x,y)  +  w(x,y) 


(B-7) 


where 


2 

V  u  =  0 

U  =  ^g(y)  ,  x  =  Ax 

U  *  ’J'yCy)  .  X  =  -Ax 

u  =  0  ,  y  =  ±Ay 

and 

2 

V  w  =  0 

w  =  0  ,  x  =  ±Ax 

w  =  »  y  =  Ay 

w  »  4»s  (x)  ,  y  =  -Ay 


(B-8) 


(B-9) 


?• 

5 


JS 


The  solution  to  u  and  w  can  simply  be  found  by  separation  of  variable  and 
added  together  to  give  the  solution  to  Problem  (1) . 


'J'l  *  I  [Cj  sinhCiyc)  +  c2  coshCp^x)]  sin  (y  -  Ay) 
n=l  ” 


+  l  [c3  sinh(umy)  +  cosh^y)]  sin  — -  (x  +  Ax) 


m=l 


where 


C,  = 


(aE  -  a:)Ei  +  tbE  -  KH  +  <4  -  4>E3 
2 Ay  sinh(pnAx) 


(B-10) 


C 


2 


<^+a:>Ei  +  <4  +  bH)E2 +  <4 +  4>e; 

2 Ay  sinh(p  Ax) 
n 


w£» to* 
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and 

ntr  _  mir 

Pn  2Ay  *  Vm  “  2Ax 

and  are  similar  to  and  c ^  but  x,  Ax,  y.  Ay,  a  and  B  are  replaced 
with  y.  Ay,  x.  Ax,  B  and  A  respectively. 

B-2.  Solution  to  Problem  (2) 

In  this  case  first  the  homogenous  equation  is  solved  and  with  the 
use  of  two  homogenous  boundary  conditions,  the  eigenvalues  are  found.  The 
solution  to  Problem  (2)  given  by  Eq.  (B-6)  is  assumed  to  be 

00 

'l'2(x.y)  =  l  cn(x)  sin  -^L  (y  +  Ay)  (B-ll) 

m=l  y 

To  find  the  unknown  function  c  (x),  Eq.  (B-ll)  is  substituted  into  the 

n 

Poisson  equation  (B-6)  resulting  in 
“  2 

l  [c"(x)  -  u  c  (x)]  sin  u  (y  +  Ay)  =  -£(x,y)  (B-12) 

m=l 

Now  the  vorticity  ?(x,y)  can  be  assumed  to  be  a  polynomial  of  x  and  y 
related  to  the  9-nodal  points  in  the  finite  differential  element.  The 
coeff  ic:?  i.-nts  of  the  polynomial  are  thus  expressed  by  the  9-nodal  values  of 
the  vorticity.  Now  that  £(x»y)  is  known  it  can  be  expanded  in  terms  of 
Fourier  sine  series  and  substitued  in  Eq.  (B-12).  On  setting  the  coeffi¬ 
cients  of  the  sine  terms  on  each  side  equal  to  each  other,  a  second  order 

ordinary  differential  equation  for  c  (x)  will  result  which  can  be  solved 

n 

accordingly.  The  homogenous  boundary  conditions  in  the  y-direction  can 
now  be  used  to  find  the  two  constants  of  integration.  The  result  of 
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Problem  (2)  thus  becomes 


ij>2  =  tc^  sinMy^x)  +  Cg  cosMy^x)  +  ax^  +  bx  +  c]  sin  (y  +  Ay) 


bAx 


C5  sinh(y  Ax) 
n 


(B-13) 


_  c  +  Ax  a 

c6  cosh(p  Ax) 
n 


a  =  -f-  (55F.'  +  a6E^  +  59E’) 
ynAy 


(B-14) 


2  (i2Ei  +  a4El  +  a8E3> 
ny 


c  =  (a  +  S3E’  +  a?E’)  +  ~  a 

y  Ay  )i 

n  7  'n 


mr 


n  2  Ay 


Note  that  a's  are  evaluated  at  the  present  vorticity  values  if  the  flow 


under  consideration  is  an  unsteady  flow. 
In  this  part  E|,  E^,  and  E^  are 


fAy 

si  =  J  A  si 

-Ay 


sin  ~  (y  +  Ay) dy 


El  = 


Ay 

’-Ay 


y  sin  ~  (y  +  Ay)dy 


■Ay  „ 

E’  =  y  sin  “  (y  +  Ay)dy 
J  '—Ay  zay 


(B-15) 
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The  FT  9-point  formula  for  the  stream  function,  i|>,  thus  is  obtained  by 
superposing  the  solutions,  and  \ given  in  Eqs.  (B-10)  er.d  (B-13) , 
and  evaluating  them  at  the  center  of  the  subregion  x  =  0,  y  *  0,  or  ^(0.0) 


if>p  =  ^l^0,0^  +  ^2^°*°^ 


I  -c2  sin  ~  +  l  c4  sin  ™ 
n=l  m=l 


+  l  (c6  +  c)  sin  -y- 
n=l 


(B-16) 


The  velocity  components  u  =  3^/ 3y  and  v  =  -3t{i/9x  are  obtained  by  the 

differentiation  of  the  stream  function  ^  and  ^  given  in  Eqs.  (B-10)  and 

(B-13)  with  respect  to  x  and  y.  Thus  the  velocity  components  at  the  center 

node  of  the  subregion  u  =  u(0.0),  v  =  v(0.0)  are 

P 


u  + 

Up  3y  3y 


c  nir  .  c  .  rn'ii 

=  l  c2  ^n  C0S  2^  +  l  C3  %  sin  T 

n=l  m=l 


r»  ,  -m  H M 

+  l  tc6  +  c]  cos  T 

n=l 


(B-17) 


!  i 


